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This report summarizes the current status and results achieved during the 
past year on research on Real-Time Flight Simulation Methodology, This work is 
a continuation of research performed during the first year on substitutional 
methods for digitization, input signal-dependent integrator approximations, and 
digital autopilot design. 

The recent Semi-Annual Report [1] was extensive and covered, in detail, 
the first six months of effort , Consequently, that report is heavily refer- 
enced herein; and only the information' necessary to provide proper perspective 
is repeated. 


I. AN INTERACTIVE SIMULATOR DESIGN PACKAGE 
FOR THE DESIGN OF REAL-TIME SIMULATORS 

1.0 Introduction 

This section describes the status of research on an interactive software 
support system which will aid the design of optimum simulation models . The 
generic type of system under study is shown in Fig. 1,0, When 'programming is 
completed, it is envisioned that the Simulator Design Package (SDP) can be used 
to evaluate a number of different standard integrator models (for example, 
Tustin, Optimum Discrete Approximation) on the basis of selectable error cri- 
teria or design an entirely new model suitable for a particular problem. In 
the latter case the model would be designed on an interactive basis, using 
selectable algorithms to find an optimal form. 

In previous work we examined a number of different substitution methods 
to determine which was most suitable under various error criteria. Based on 
these results, a number of substitution formulas have been chosen for inclusion 
in SDP , Consequently , most of the work during the first six months concen- 
trated on the evaluation of optimization algorithms and discrete representa- 
tions, The second six months have been used to develop a substantial amount of 
the software framework of SDP. This includes subroutines for iterative designs 
of simulation models as well as a rudimentary graphics package. 

' This section contains three sub-sections which report on different facets 
of the development of SDP , In the first sub-section we are concerned with 
assuming a discrete representation for a given continuous transfer function and 
then iterating to solve for optimal values of the parameters concerned. Pre- 
liminary results of time domain optimization are also discussed. In the second 
subrsection a similar effort is reported in which a form is assumed for an 
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integration operator, and then a random search method is used to determine the 
optimum parameters. In the third sub-section examples of graphical results 
obtained with SDP are presented, 

1.1.1 Frequency Domain Optimization 

As indicated in previous reports, attempts have been made to obtain a 
digital simulation model in a way that minimizes the frequency^domain error. 

We have discussed the general form of the digital system, the error criteria^ 
and the gradient technique used. Implicit formuli for gradients were derived 
so that programs can be written to evaluate gradients necessary for the numeri- 
cal technique being used. 

Referring to the 1976 Semi-Annual Report [1], we have the following 
results : 

Error: 

, M 

E = S H(joj^) - ^ (1-1) 

' ^ =■ 1 , ■' ■ 

Form of digital system; 

n . (1 + (1 + 

H(z) = A jjQp /2 ~~~ nF ~ (1-2) 

n (1 + 2b_cos2b, z-1 . ,2 ^-2, n (1 + b, z 
k = l \-l \ = NCP+1 

where NZ = number of zeros 

NP = number of poles 
NCP = number of complex poles 
k E largest integer NZ/2 

“ 0 if NZ = 2k (even number of zeros) 
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0 if NZ - 2k + l(odd number of zeros) 

We also discussed the constraints necessary to limit the time, domain 
error and to ensure stability. 

Three FORTRAN programs have been written to accomplish the optimization 
in which the Fletcher-Powell method [2] was used. Since the algorithm is not 
critically important and the Fletcher-Powell method has proved to be success- 
ful, no other numerical methods were investigated. The process of obtaining 
the optimal digital form will be implemented off-line; so the speed-versus- 
accuracy criterion is no longer important in choosing the numerical algorithm. 
The Fletcher-Powell method converges rapidly in this particular application. 

Application of the Frequency Optimization (FO) Method 

Only recently were we able to use our programs to perform some prelimi- 
nary comparisons. This is due partly to the complexity involved in program- 
ming, The programs were written in such a fashion that they can be incorpor- 
ated into the SDP in the future. Only first-order systems have been investi- 
gated, to date, Higher-order systems will be considered in future reports. 

The first-order system under study is: 

a-3) 

The selected sampling interval is T =11/10 sees. This is rather large, 
but it will "provide a good comparison with the Tustin method for large sampling 
intervals. Most existing simulation methods work fairly well with small sampl- 
ing periods but suffer severely at larger sampling periods. 

The frequency response is given by 


4 




H(j(o) = — (1-4) 
1 + u» 

The actual response of this system will be compared with the responses of the 
simulation system for a step, ramp, and sinusoidal inputs. 

Several different optimized systems were obtained and their performances 
compared to the above system. A distinct characteristic of this optimization 
procedure is that the user can arbitrarily specify the order of the simulation 
system, i.e, , one can obtain a first-, second-, or third-order digital model 
for the simulated first-order system. As results will show, it is usually 
better to use the digital model of the same, or one order higher, than the Con- 
tinuous system. Furthermore, if the continuous system’s poles are all real, 
then the user should choose his digital model to have all real poles in the z- 
plane . This is intuitively obvious, since all real poles in the s-plane map 
into real poles in the z-plane. 

As discussed in the semi-annual report, in order to limit the time domain 
error we must place constraints on the pole locations of the digital system in 
the following fashion: 

z-plane poles = exp [T (s-plane poles)] (1-5) 

or in a less restricted form: 

I I I sTi T (Real s) ,, 

|z| = |e 1 = e ^ , (1-6) 

This form determines the radius of a circle (in the z-plane) . All poles (both 
real and complex) must be constrained to lie within this circle. We also allow 
this radius to approach unity in order to show that the constraint does, 
indeed , reduce time-domain error . 


The following digital systems are models of the first-order continuous 
system which we simulate: 

(A) Poles of digital systems are restricted to lie within the 

-T 

circle of radius e = .730403 


1 -Order Model: 


H(z) - K 


1 + .9866292 
1 - .713722z' 


where K = 1.144102 


2 -Order Model: 


H(z) = K 


1 + 2.42315z~^ + .518598z~^ 

(1 - .729328z‘^)(l + .730403z“^> 


where K = .118823 

(B) Poles of digital systems are restricted to lie within the 
circle of radius .9 


1 -Order Model: 


H(z) = K 


1 4- .9866292 
1 - .713722z’ 


where K =• , 144102 


2^'^'T-Order Model: 


1 + 2.78393z"^ + .7439267,"^ 


728308Z ) (1 + .90000Z ) 



where K = .114034 


(C) Poles of digital systems are restricted to lie within the 
circle of radius . 999 

l®*^-0rder Model: 


H(z) = K 


1 + V986629Z 


-1 


1 - .713722Z 


-1 


where K = .144102 


2^‘^-Order Model ; 


H(z) = K 


(1 + 3.01207z~^ -H .885306z~^) 
(1 - .727803z"^)(l + .999z"^) 


( 1 - 11 ) 


( 1 - 12 ) 


where K = .111105 

All K are chosen so that the steady-state gain is 1, which is the steady-state 
gain of the continuous system. 

Notice that all first-order models are the same, regardless of the con- 
straint, This is possible, because we are simulating a first-order system. 

For second-border models the constraint plays a more important role. The opti- 
mization algorithm tends to place poles as far as possible within the con- 
straint, This results in poles at -,9 (for | Poles] ^ ,9 constraint) and at 
t, 999 (for [Poles] ^ ,999 constraint) . These poles dominate the response of 
the system and^ as a result, give more error, 

The following figures show some of the advantages and disadvantages of 
the FO method. The method is far better than Tustin, especially at high sampl- 
ing periods. As stated before, the digital systems are designed for a sampling 
period of T = 11/10 secs. This is rather large for the Tustin method, whose 
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Figure 1-1 Comparison of Exact, Tustln, and 1st Order FO Responses 



Figure 1-2 Comparison of Exact, Tustin, and 1st Order FO 
Responses (Constraint Relaxed) 




Figure 1-4 Comparison of Exact, Tustln, and FO Responses Ramp Input 




performance is best at smaller sampling periods. As we might expect, relaxing 
the pole magnitude constraint results in slight oscillations in the simulated 
responses. Empirical results also show that better time-domain accuracy can be 
obtained by using a digital system of one order higher than the simulated con- 
tinuous system. This may result in slightly more computation time, but it is 
not particularly critical. The digitaT system can be designed with high sampl- 
ing periods to accomodate the real-time constraint. 

One disadvantage of the FO method is that it is not very flexible in 
terms of sampling period. Once a digital system has been obtained, the sampl- 
ing period cannot be changed without degrading the simulation performance. If 
the designer desires a model for a different sampling period, he has to start 
the design process again. This characteristic is shown in Fig. 1-6, where the 
digital system had been designed f or T = II/IO secs . ; but the simulation was 
performed at T =?= ,1 seC, For sampling periods other than the one that has been 
designed for, the time-domain error increases substantially . Fortunately, this 
drawback is not particularly severe. The design procedure depends almost com- 
pletely on the digital computer, and there is little calculation to be per- 
formed by the designer. 

All in all, the FO method seems very promising from preliminary results 
obtained to date , Only first-order systems have been tried , and no conclusive 
results can be derived from this limited empirical data. The next task will 
include higher^order systems with a variety of inputs. 

1.1.2 Time Domain Optimization 

The preceding discussion was centered around the time-domain response 
resulting from a frequency domain optimi zation , Although there is a unique 
relationship between time and frequency response, minimizing frequency-domain 
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Figure l-"6 Goinparison of Exact, FO, and Tus tin Responses for T = .Isec 


error does not necessarily make time-domain error a minimum. Our ultimate goal 
is to find an optimal time-domain simulation. 

The state-of-the-art of optimal time-domain simulation is virtually non- 
existent. The biggest obstacle seems to be the definition of a general per^ 
formance index whose value must be tninimizad to yield a desired digital system. 
Unlike the frequency domain performaixce index^ the time-domain error is input 
dependent. This feature makes it difficult, If not impossible^ to obtain a 
general performance index that will suit all inputs which vary from application 
to application. Presented below are some of the results published in the lit- 
erature. Most of the methods presented are far from being applicable to our 
purposes. These approaches, however, could be used for our application, since 
the performance index depends only on the numerical values at the sampling 
instant of the ideal response. 

J. A, Athanassopoulos and A. D. Warren [3] proposed the following form of 
the digital system: 



m 


E a z 
i = 0 ■ 

H(a, z ) - 


-1 


n 

E z 
i = 0 ^ 


-1 


(1-13) 


where 


a 








(1-14) 


The time-domain response of the above system is a function of the coefficient 
vector a whose determination depends on the performance index . 

Now, let be a specific value for the response at 

t = tj^(k = 1, 2, • • • ), let s be a set of all real a such that poles of 
H(a, z b lie outside the unit circle, and let R(a; k) be the response of 
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interest. In this application R(a, k) is chosen to be the time domain response 
y(a> tjj.) = 

The design problem is then of the following min-max form: 

Given the specification set =‘1, • • * } and the maximum 

allowable degree of H (i.e. , maximum numerator and denominator 
order), find the coefficient vector a = a* such that: 

' max|R(a*, k) - r^l ^ max|R(a, k) - tj^l (1-15) 

k > k 

Notice that a*e s in order that the digital system be stable. 

The above problem is converted to the form of a mathematical program by 

using an additional variable n. Thus, in order to solve (1-15) a vector 
■ T 

a - [ot. n] which minimizes n subject to the constraints 


4'^(o, k) = R(a, k) + ® 


(1-16) 


^ 2 ( 01 . k) = -R(a, k) + Y2kH + ^ 0 


(1-17) 


must be determined. This can be easily seen in Fig. 1-7. 

The goal is to minimize the width ® zone subject to the 

constraint that both R(a, k) and rj^ lie in the zone for all k. Various meth- 
ods, including a penalty function, can be used if the constraints are highly 
nonlinear functions of the variables, 

More studies must be directed to this method before it can be directly 
applied to our application. The constraint is only that poles must be inside 
the unit circle. This is not enough to obtain a good representation of the 
continuous system, since the optimization procedure will tend to place poles at 
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farthest possible location. 

Other approaches proposed by Burrus and Parks [4] , and 
Steiglitz [5] have appeared in the literature. Unfortunately, no stability 
considerations were given for these methods; and, as a result, they are not 
practical for our research. With these discussions in mind, we propose the 
following approach to our research in tlwe-domain optimization. 

The digital system has the same form as used in the FO method, i.e.. 


Y(z) 

U(z) 


= H(z) = A 


. k = 1 : : ' V ^ 

MPP / 9 NP 

n (1+ 2b2.cosb2, , z"^ + n (1 + b z"h 

k = 1 ^ k = NCP+1 


(1-18) 


The advantages of this digital form were discussed in the semi-annual report. 
The error function is defined in time domain as: 

N - - ■ 

E = S |Y - Y(nT)| (1-19) 

n n . . . ■ ■ . . ■ ■ 

. n = 0 

where Y is the exact response and Y(nT) is the simulation response from (1-18) 

V „n. , / .. 

for a given input. This error function is input dependent; and, in order to 
achieve maximum generality, we must divide inputs into different classes and 
design the digital system accordingly. In other words, a prototype digital 
system is designed for each class of input; and the designer must choose 
between several systems for the best one under a given circumstance. This 
approach is an extension of the Sage ' s Optimum Discrete Approximation which 
only uses the step and ramp inputs as test inputs . We can extend this further 
and use a variety of inputs with the aid of the computer and the optimization 
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algorithm that Is already available. 

1»2 Digital Simulation and Optimization via Random Search Techniques 
1»2.1 Introduction 

This section continues discussion on our work involving techniques for 
the development of a discrete time integration operator to be used in the simu- 
lation process. The integration operator can be optimized for a particular 
system subjected to a set of specified Inputs. The class of systems being 

investigated are those which can be represented by a set of state equations. A 

■> 

discrete time integration operator with certain free parameters is hypothe- 
sized. An adaptive random search optimization (ARSO) technique is used to find 
the .optimum values for these parameters, Examples are presented to show the 
effectiveness of this technique. 

1.2.2 Integration Operator 

The class of systems being investigated are those which can be repre- 
sented by the set of state equations 

X = f(x, u) (1-20) 

where x is the n x 1 state vector, u is for the r x 1 control vector, and f is 
the set of n functions, typically nonlinear. 

Figure l-8a is a block diagram of the mathematical relations in Eqn. 
(1^20) , The vectors x and u are acted upon by the functional relations 
f(x, u) , producing the vector x. Figure l-8b is a block diagram of a discrete 
approximation to the continuous time system. The control vector u is assumed 
to be sampled at a uniform rate, producing the input samples u(k). The 
equations 
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Figure 1-8 Digital Simulation of Nonlinear Systems 
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x(k) = f[x(k), u(k)] 


( 1 - 21 ) 


are in the same form as those representing the continuous time system. A 
discrete-time integration operator of the following form is chosen. 


F(z) 


N 

T Z X.z^ 
j = 0 ^ 


N 


z 


(z - 1) 


N 

T(Xq + X^z + • • • + AN^ ) 

__ N + 1 N ^ 

z - z 


F(z"^) 




- 1 


Z-2 + 


j. ^ -N . , -(N + 1), 

+ A^z + AqZ ^ '] 


1 - z 


-1 


( 1 - 22 ) 


where T is the sampling period, and the A 's are a set of free parameters, the 
values of which are to be optimized. This operator yields a realizable simula- 
tion, since the power of the denominator is always one greater than that of the 
numerator. The pole at z = 1 corresponds to a pole at the origin in the com- 
plex s-plane , and the N*’^-order pole at the origin in the z-plane corresponds 
to an N^^-order pole at negative infinity in the s-plane [; 6] . Therefore, the 
transient response of the poles added at z = 0 to make the operator closed-loop 
realizable will decay quickly. The state equations are now of the form 

x(k) = f[x(k), u(k)] 


and 

X(k + 1) = x(k) + T[AjjX(k) + Ajj_^x(k - 1) + • • • + AgX(k - N)] (1-23) 

Equation (1-23) can be thought of as a polynomial approximation to the value of 
a function at point (k + 1) based on its value at point (k) and the value of 
its derivative at point (k) and preceding. 
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The free parameters in F(z) are optimized using an idealized model form 
of a model reference adaptive control system [7]. Figure 1-9 is a block dia- 
gram of this configuration* 

1.2.3 Optimization Technique 

The perturbation of parameters in F(z) is controlled by an Adaptive 
Random Search Optimization technique (ARSO) which was described in our previous 
report [1] . 

The unknown parameters are perturbed in the following manner: 

+ 1) “ A* + (Sx^a + 1) ^ (1-24) 

th ' 

where A(j + 1) is the new value of the i parameter, X* is the "best-to-date" 
til 

Value of the 1' parameter; that is, the value. of X^ when the minlmum-to-date 
value of the J vector was calculated and dX^^Cj + 1) is the random perturbation 

t*'h 

for the i parameter. This. is equivalent to the perturbation scheme shown 
belo% 

' x^(j + 1) = X^(J) - a(j)5X^(j) + 6X^(j +1) 

where . 

a(j) - 0, I* 

-1. if Jy, -1) 

j is the performance vector in the trial, and is the smallest 

performance vector ohtained through (j - 1) trials. The coefficient a(j) is 
used to negate the effect of an unsuccessful trial. The perturbation is 
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calculated as follows : 


<SAi(j) = V^a) +/3a^(j)[2RND(0) - 1] (1-25) 

th 2 

where y^(j) is the current value for the mean of the i random variable, a (j) 

is the current value for the variance, and RND(O) is a uniformly distributed 

random variable on the interval [0, 1]. Equation (1-25) produces a random 

number from a unifomly distributed random variable with mean y and variance 
2 

a • Stability considerations may place constraints on the parameter calculated 
in (1-24). If a particular perturbation places a value outside its limit 
for stability, the value may be moved deterministically inside the limit, or 
another random perturbation may be tried. In simulating complex systems, how- 
ever, it may be difficult to a priori determine the stability limits for the 
coefficients. In this case one or more of the state variables may be monitored 
during the simulation; and, if they exceed reasonable values, the trial may be 
aborted. This saves computation time and may prevent the entire program from 
being terminated due to overflow. 

When a particular trial is successful, that is, 

Ji 1 J| 1 £ i 1 n 

where J* is the minimum value of the i element of J, the means and variances 
of the random variables are updated. The mean is calculated as follows : 

- m^ (1-26) 

where 
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c . (1) c . (2) c (3) C (4) G (5) 
m. = -i— + + -f— + 

X 2 4 8 16 16 


(1-27) 


The c^*s are past values of A*; that iSj previous values of the "best-to-date*^ 

parameters. The most recent value of A* is c.(l). Since the most recent best 

i X 

value corresponds to a smaller index than a previous best value, (1-26) 
tends to select the most favorable direction for the next perturbation. 

The variance for the distribution is determined using the following argu- 
ment. The perturbation for each parameter is a random variable with a uniform 

■ ■ ■ ■ '■2 ' ' 

probability density function with mean y and variance a . Relating these 

moments to the end po^^hts of the function (a, b) yields 


y = (b + a)/ 2 


(l-28a) 


= (b - a)^/12 


(l-28b) 


If the mean and one end point is known, the other end point and variance can be 
calculated. That is, if y and a are knov/n, then 


b = 2y - a 


(1-2 9 a) 


= (y - a)^/3 


(l-29b) 


and if y and b are known, then 


a - 2y - a 


(l-30a) 


a 2 = (b - y)-/3 


(l-30b) 


Figure 1-10 illustrates two possible situations with a scalar cost fi 


tion and a single parameter after five successes, In the top figure, c (2) 
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should be the end point (a) to ensure that the optimum parameter value lies 
within the search area. This gives 

a = c(2) - c(l) (1-31) 

where c(l) , the current best value, is taken as the zero point on the density 
function. Substituting this into (1-28) yields 

0^ = [y - c(2) + c(l)]^/3 (1-32) 

In the bottom figure c(2) should be the upper end point (b) . Thus, 

b = c(2) - c(l) (1-33) 

and 

0^ = [c(2) - c(l) - 

Equations (1-32) and (1-34) yield the same numerical value for the variance, 
since the sign difference is lost when the numerator term is squared. There- 
fore, (1-34) , combined with the mean from (1-26) , provides the necessary 
data for computing the perturbations. IThen a large number of consecutive 
failures are generated, the means and variances are set to deterministic 
values. This allows the search technique to more fully explore the parameter 
space, about which little is known beforehand. If no improvements are obtained 
after a set number of failures , the search is terminated or restarted with a 
different set of initial conditions. 

1.2.4 Results 

One of the problems considered is the simulation of the nonlinear equa- 
tions of an aircraft. These equations are shown below: 
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X = V cos (y) 


. , . d 2 ^ T , . 

V = -g sin(Y) “ ^ ^ cos (a) 

Y = (1/v) (-g cos(y) ^ ^ sin (a)) (1-35) 

' ■ ■ a =■ 0 ) , ' . . . . ■ . . 

0 ) = 1.311 u ~ .806 0 ) - 1.311 a 

xvhere x is the horizontal displaGeiuent , h the vertical displacement, v the 
total velocity, y the flight path angle, a the angle of attack, u the control 
elevator deflection, and o) the rate of change of the angle of attack. 

Each of the si^x states generates an element of the performance index. 

The mean squared error between the approximate value and a value obtained by 
Runge>-Kutta integration is used, The elements are; 

J1 = range error (feet) 

J2 - altitude error (feet) 

J3 = velocity error (feet/sec) 

J4 ~ error in flight path angle (radians) 

J5 = error in angle of attack (radians) 

J6 = error in angle of attack rate(rad/sec) 

The RMS results of ARSO optimization using one, two, and three parameters 
are shown in Table 1~1. The sample period for each of these is 0 ,5 seconds , 

Also shown in the table are the results using the forward difference operator 
and the Milne-Reynolds predict or‘-corrector method. The Milne-Reynolds results 
are for a sampling period of 0.25 second , as both it and the delayed Tus tin 
are unstable for T = 0 , 5 second . The total amount of computation per second 
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of simulation time was four times greater for Milne-Reynolds than for the dis- 
crete operators obtained via ARSO. 


Table 1.1 


Preliminary ARSO Results 


ARSO-1 ARSO-2 ARSO-3 


J1 


218 


6.796 

2.121 

J2 


4.39 


3.63 

3.385 

J3 


.305 


.216 

.205 

84 


8.76.10~‘^ 


4.31.10“^ 

4.16.10"^ 

85 


3.0.10“^ 


1.39.10"^ 1.37.10 ^ 

J6 


3.4.10“^ 


1.59.10"^ 

1.57. 10"^ 


■ 

Fwd. Diff. 

Milne Reynolds 


J1 

3,69 


5.46 X 10 ^ 


82 

; 6.09 


I.’SI ,x 10”^ 


83 

.371 


■ 3.22 X 10~^ 


84 

7.4 X 

1 

0 

4.6 X 10"^^ 


85 

: 3.1 X 

10"^ 

2.7 X lO”^ 


86 

3 . 5 X 

10“^ 

■ 1.4 X 10~^ 


ARSO-l: = 1 . Q127_ 

ARSO-2; = 1.20793; A2= “•208282 

ARSO-3 : >= 1. 20999 • A^ = - . 200159 ; A 3 = -9 . 83883 .10“^ 


The aircraft maneuver for the data in Table 1.1 was a climb followed by a 
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dive with a total range of approximately 5, 6 miles. To evaluate the sensitiv- 
ity of the cost function to changes in input, a second maneuver, a dive, was 
run. The total range was about the same as the first run, The results of this 
maneuver for ARSO-2 and ARSO-3 are shown in Table 1.2. 


Table 1.2 


ARSO-2 


ARSO-3 


J1 

8.96 

2.82 

J2 

2.48 

2.38 

J3 

.117 

. 113 

J4 

4.13.10"^ 

3.91.10"'^ 

J5 

9.49.10"^ 

9.38.10“^ 

J6 

1.1. 10~^ 

1.09,10"^ 


Using the ARSO-2 and ARSO-3 coefficients from maneuver one as initial 
conditions, the parameters were then optimized for input two. Table 1,3 pre^ 
sents the results of this optimization, and Table 1.4 lists the cost functions 
when maneuver one is run with the coefficients optimized for maneuver two. 
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Table 1.3 



ARSO-2 

ARSO-3 

J1 

1.32 

.653 

J2 

2.48 

2.35 

J3 

.117 

.111 

J4 

4.04. 10”'^ 

3.91,10‘ 

J5 

9. 51 . 10 “^ 

9.3.10"' 

J6 

1.1. lO'^ 

1.08.10' 



Table 1.4 



ARSO-2 

ARSO-3 

J1 

3.99 

2.92 

J2 

■ : 3.53 - 

3.32 

J3 

.212 

.200 

J4 

:4> 34 . 10 "^ 

4.11. lO' 

J5 

1.39. 10~^ 

1,35.10' 

J6 

1 . 60 . 10 “^ 

1 , 55 . 10 ' 


ARSO-2: =1.20807; ^2 = -0.207852 

ARSO-3: = 1.21322; A 2 " -•201530; A 3 = -1.15446.10 


As can be seen from the tables, increasing the number of parameters 
decreases each element of the cost function in each case, This is because the 
higher-order polynomial approximation of (1-24) can better represent the 
actual function. It should be noted that, although the coefficients are input 
dependent, neither they nor the cost function change significantly when differ 
ent maneuvers are executed. As an aid in evaluating the results of Tables 1,1 
1.4, Table 1.5 lists the cost function for ARSO-2 from Table 1.1 as a percent- 
age of the dynamic range for the respective state variables. 


Table 1.5 



ARSO-2 

Dynamic Range 

Percei 

J1 

J2 

6,796 

3.63 

29,772 

0,02: 

1,19 

J3 

.216 

—4 

14.318 

1,51 


4,31.10 

^2 

.0525 

0,82 

J5 

1,39,10 

1.0813 

1,29 

J6 

1.59,10 

,2607 

6,1 


Since the performance vector does not allow an increase in value for any of the 
elements, the ARSO technique is dependent on initial parameter values. Thus, 
several different sets may be used before the final pararaeter values are ^ 
determined. 
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1.2.5 Further Study 

It is desired to apply the ARSO technique to a twelth-order set of equa- 
tions representing a six degree of freedom aircraft. This will be a signifi- 
cant step forward in determining the general applicability of this optimization 
technique . 

1.3 Some Options Available SDF 

The existing programs of the SDF allow the user to siniulate a continuous 
system up to 10 order using the "Tustin and Optimum Discrete Approximation 
methods . Results are compared with the fourth-order Runge*-Kutta method to 
determine ' the error. The Runge-Kutta method is applied at a smaller sampling 
period to achieve ’-ideal” response. 

Various plots can now be obtained; For comparison purposes ■ the time 
response for each method is plotted against the ideal response for a given sam^ 
pling period or each method is plotted for different sampling periods. Regard- 
less of the type of plot the user may desire, there always is an ’’ideal’f plot 
from which he can observe how well a given siraulation method performs , 

At the end of each simulation run, the user can also obtain an error 
plot . This error is based on the deviation from the Runge-Kutta method and is 
plotted against each sampling period that the user has previously specified. 

To illustrate all the options available, two linear continuous systems 
were simulated. In all cases part (a) of a figure illustrates man/machine 
interactions , while part (b) shows the graphical results. The first example is 
a seventh-order autopilot transfer function. Step and sinusoidal inputs result 
in plots in Fig. 1-11 and Fig. 1-12 with sampling period of 
T = ,15708(40 rad/sec) . The second example is a fourth-order continuous 
system. Figure 1-13 shows plots of responses obtained from the Tus tin. Sage, 
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Figure 1-11 Step Response Using T - .15708 
Illustrating Tustin and Sage 
Methods for 


F(S) 


36 S^ + 2.725 + 2.31 S +1.65 S^ +7.255+81 

. — • — — ... ■ ■■ ■ — • — — — • ^ ■■ • . ' . • I., ,1 I — 

S^ + 8.45 +36 S^ + 5.62 S + 3,1 S + .62 1.125S^ + 13.33 S +81 
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INPUT NUrJERATOR AMD DENOrtlNATOR ORDERS 

5>7:. : : ' 

muI’^erator coefficients in ascending order 

17970.9, 1 1933. a;37l9. 74,358. 78, 33, 
rESOfllNATOP CCEFFICENTS IN ASCENDING ORDER 
-i9S:t. 83 ,19049. 1,23383. 3, 10433. 1,34. 6 1, 6, 340. 468,26. 4889,1. 
INPUT INITIAL AND FINAL TINES 
0> 15,708 

HOU NANV SAMPLE PERIODS'? 

I ^ 

SANPLE PERI0D5-GREATER THAN.314E-01 
.15708 

RESPONSE LISTED? YE$*0,NO*1 
i 

RESPONSE P1.0TTED? YES-0,NO-1 
0 

PERFORMANCE INDEX; 0-f1SE,l-MAE,3-riTAE 
0 

INPUT E 0-STEP, I -RAMP, 2-SlNE 
0 

AMPLITUDE 


.2493E407 

.ed00£40« 

.0000E+00 

.0000Et00 

SAGE DENOMINATOR COEFFICIENTS 

-.434E4^06 

•4730E+07 

-.230E+08 

.6313E^08 

-.105E>e9 

.1042E+89 

-.576E'*-e8 

.1356E408 

.0000E+00 

i- .0d00£4.00 

^ .ee00E+00 


P.I. SAMP. PER. 

.1344E-02 .157E+ea 


TUSTIN NUMERATOR COEFFICIENTS 
-.734E+07 
. 1975Efe8 
-.158E+08 
-.I26E+08 
.5067E^08 
-.364E>98 
-.268£^08 
.2999E+08 
.e«00E>00 
.00006^60 
.0000E+00 


TYRE EGO FOR RESPONSE PLOT 

sinuL ; PAUSE eee? 
sinuL SUSP 




TUSTIN DENOMINATOR COEFFICIENTS 

-.861E+07 

,8ll8E+08 

-.36SE>09 

.1003£*10 

-.178E+10 

.2016E+10 

-.133E+10 

.3713E409 

.0090E+00 

.0000Ef00 

.0000E+00 


P.I. SAMP. PER. 

.4868E-03 .157E4>00 


SAGE NUMERATOR COEFFICIENTS 

.0000E400 

.0e00E400 

-.335E+06 

.2262E+07 

-.666E+07 

.1O34E+08 

-.799E+07 


Figure 1-11 (a) 


3 . 5 ^ 



Eigure 1~ 



F(S) = 


Figure 1-12 Sinusoidal Response (w = 1) for T= .15708 
Illustrating Tustin and Sage Methods for 

36 S^+ 2.72 S + 2.31 S + 1.65 s2+7.25S+ 

— M i i« - - ........ • I., • — ■ ' ■ — .■■■■ ■ ■ 

S^ > 8.4S + 36 s2 + 5.62 S + 3.1 S + .62 1.125S^ + 13.33 


81 

+ 81 


INPUT HUrtERATOP AND I>€NOniNATOR ORDERS 
5.7 

HCMERaTOR coefficients in ascending ORDER 
1 1S32. 9. 17970. 9. 11923. 2/371S. 74. 358. 72.33- 
rEnOfllNATOR C0EFFICENT5 IN ASCENDING ORDER 
4551.82.19049.1.23363.8. 10433. 1,2461. 6. 3*46.468. 26. 
1NP.:T initial and FINAL TINES 
0. 15.708 

HOU NANV SANPLE PERIODS? 

X 

SAMPLE PERIODS-CREATER THAH.314E-01 
.15708 

RESPONSE LISTED^ YES-0.NO-1 

1 /, . . 

RESPONSE PLOTTED? YES-0.NO’»1 
0 

FERrORMANCE :NDEXi0-f1SE.l-nAE.2-nTA£ 

0 

INPUT: 0-STEP. 1-RANP. 2-SlNE 

2 

ANPLITUDE AND FREQ. -RAD/SEC 

1.1 


TUSTIN NUMERATOR COEFFICIENTS 

-.724E407 

.1975E^eS 

-.158E>08 

“.lEGE+es 

.5067E408 

-.364E408 

-.268E^08 

.2999E>08 

•0e00E400 

.0000E^00 

-00e0£+e0 

TUSTIN DENOMINATOR COEFFICIENTS 

-.861E+07 

.8118E408 

-.365E409 

.1003E+10 

-.178E410 

.20166410 

-.132E+10 

,.3713E+eS 

.0e00£4ee 

.08006400 

.00eeE400 


P.I. SAMP. PER. 

.7834E-03 .157E400 


SAGE NUMERATOR COEFFICIENTS 

.0000E400 

.00e0e400 

-.335E+0S 

.2H62E407 

-.666E407 

.i024E408 

”.799E+07 


.2493E407 

.0000E400 

. 0000 E *»00 

.00006400 

,1. SAGE DENOMINATOR COEFFICIENTS 

-.4246406 
.4730E407 
-.2306408 
.6^13E*^e8 
-.105E409 
.1042E409 
-.576E468 
.1356E408 
.00006^00 
.0000E400 
.0000E400 


P.I. SANP. PER. 

.44e9E-03 .157E400 


TYPE :G0 FOR RESPONSE PLOT 
SIHUL : PAUSE 0007 
SIMUL SUSP 



Eigure l-12(a) 


EXACT -SOLID LINE 
TUSTIH • XXX 
SAGE • +++ 


figure 1-12 (b) 


Figure 1-13 Step Response for T = .05 Illustrating 
Tustin and Sage Methods for 


F(S) 


S^ + 5.5S + 2.5 

+ 5s3 + 9.5S^ + 8S + 2.5 



INPUT HUNER*^TOR PHD DENOMIHATOR ORDERS 
2.4 

NUniRATOR COEFFICIENTS IN ASCENDING ORDER 
2.5;5.5/l 

lENOniNpTOR C0EFFICENT5 IN ASCENDING ORDER 

INPUT INITIAL AND FINAL TINES 

0,10 

HOU NANV SAMPLE PERIODS? 

1 

SAMPLE PERIODS-GREATER THAN-20OE-01 
.05 

RESPONSE LISTED? VES-0.NO‘l 
I 

RESPONSE PLOTTED? VES*0,NO-1 

e 

PERFORnANCE INDEX : 0-NSE, l-nAE,2-NTAE 

e 

input: 0-STEP, 1-RANP, 2-SINE 
0 

ANPLITUDE 

T 


.eoeoE^eo 

• 0000E4>00 

•eooeE^oe 

.eoo0E4>00 

SAGE DEHOniHATOR COEFFICIENTS 

.1600E+C6 

-.G80E+06 

.IO84E+07 

-.768E+66 

.2O40E+06 

•0000E+e0 

.0000E+e0 

.eoooE^ee 

.00O0E+00 

.ooeeE't-oe 


P.I. SANP. PER. 

.1327E-03 .50,0E-01 


TUSTIH NUNERATOR COEFFICIENTS 
.l38^E-^e4 
-.430E+03 
-.318ET04 
.45eeE>03 
• 1822E+04 
.ee00E^00 
.0e00E^00 
4S .0000E>06 

H .0eeeE+e0 

.0000E+0e 

.eeeoE+ee 


TVPE :G0 FOR RESPONSE PLOT 
SinUL : PAUSE 0007 
SINUL SUSP 



TUSTIN DENOfllNATOR COEFFICIENTS 

,2255E+07 

-.96ee+07 

.1533E+08 

-.lesE+es 

.2896E+07 

.000OE+00 

.000e£-^00 

.ee0eE+00 

.0000E+00 

.0e00E+00 

.ee00€+00 


P.I. 


SANP. PER. 


.I262E-03 


.seoE-di 


SAGE HUHERATOR COEFriCIENTS 

.0000E^00 

•00C0e+00 

.4000£^03 

-.910E+03 

.5125E+03 

.00eeE+00 

,0000E+00 


Figure 1-13 (a) 



Figure 1-13 (b) 


and fourth-order Runge-Kutta methods for a sampling period of T = .05 sec. ' 
Figures 1-14 through 1-18 are plots of responses of a given method (in this 
case, Tustin) with different sampling periods. The mean-square-error for each 
sampling period is calculated and is plotted in Fig. 1-19, 



43 



Figures 1^14 - 1-18 


Step Response F or 
T = .2, .4, .6, .8 
and 1.0, Respectively 
Illustrating Tustin 
Method for 


4- 5.5 S + 2.5 

4- 5 S^ + 9 . 58^ 4- 8 S 4- 2 . 5 


F(S) 


IHPirr NUhERATOR AND DEHOfllHATOR ORCERS 
2,4 ■: 

riUMERATOR COErriCIENTS IN ASCENDING ORDER 
£.5/5.5. 1 

DEiSOniNATOR COEFFICENTS IN ASCENDING ORDER 
2.S.S.9.5.5.1 

INFUT INITIAL AND FINAL TIMES 
0.,i0. 

Hou MANV Sample periods'? 

5 

SAMPLE PERI0D5-CREATER THAH.20OE-ei 
.2. .4.. 6.. a. I . 

RESPONSE LISTED? VES-O.NO-l 
1 

RESPONSE PLOTTED^ VES-0,NO-1 
0 

P£Rfi:)RMAHCE index:: 0-flSE, 1-MAE/2-MTAE 
0 

METHOD: 0-TUSTIN; 1-SAGE; 2-IBM 
0 

INPUT: 0-STEP, l-RAf1P, 2-SlNE 
0 

AMPLITUDE 

1 ; - 


TUSTIN NUMERATOR COEFFICIENTS 
.4750E^^02 
-1.00€^'02 
L85E+03 
.1200E^03 
.1575E*^03 
.0000E+00 

.0000€+0a 

.^0E+00 

.0000C+00 

,0000£+00 

-0000E+00 

TUSTIN DENOMINATOR COEFFICIENTS 

.5872E+04 

-.301E+05 

.5811E+05 

-.498E>05 

.1603E+05 

.0000E+00 

.0000E+00 

.0000E+00 

.0000E^00 

.0000E+00 

•0000E400 


TYPE ICO FOR RESPONSE PLOT 
SIMUL 1 pause 0003 
SlfHJL SUSP 





Figure l-14(a) 


0 ^ 


a.5«- 

H. 25’ 

a.oe- 

I . 75 ' 
1 . 50 ” 
1 . 25 ” 
1 . 00 ” 

.75 ■ 
. 50 ” 
. 25 ” 
. 00 ” 


SrtnpLE PERIOD •.aeeaE+ee 


EXACT • SOLID LIME 
TUSTIM ■ XXX 
SAGE • +++ 


TYPE IGO TO CpMTINUE 
SIMUL J P AUBE 0004 

simjL SUSP f 


Figure 1-14 (b) 


TUSTIN NUMERATOR COEFFICIEMTS 

.i»O00E>e0 

-.-»50E+0a 

-.350E>d2 

,65O0E+02 

.5500E>d2 

•0OO0E+00 

.OOO0E+0d 

.0oooE+eo 

.o0ooE>ee 

v00O0E>oe 

^oooeE^ee 

TUSTIN DENOrtiHATOR COEFFICIENTS 
-a0O0£>03 
>.132ET04 
.3290£^O4 
-•366E+e4 
.1S30E404 
.0000E+00 
.0oe0E^00 

.0000E^00 

.0000E>00 

.0000E400 

.0000E^“00 


TVPE tCO FOR RESPONSE PLOT 
SinuL : PAUSE 0003 
SinUL SUSP 


• 


Figure 1-15 (a) 



X 



EXACT • SOLID tlH€ 


SAtiPLE PERIOD -.^eeaE+e® tustin • xxx 

SAGE ■ ♦♦♦ 

TYPE SCO TO CpHTlHUE 
SinUL : PAUSE 
sincL SUSP 




Figure l-15(b) 


TU5T1N HUrtERATOR CbEFriCIENTS 

-.267E4ea 

-.722E+01 

.4667E+e2 

.319-4E^^02 

.0OO0E40O 

.00O0E406 

.OO0OE^00 

.O000E400 

.eeoet^ee 

.O000E+00 

TUSTIN DENOniNATOR COEFFrCIENTS 

.1966E+02 

-.167E^03 

. 54-46E+03 

-.301E>03 

.^134E>03 

.00O0E'fO0 . 

.0000E'»’00 

.d000E+00 

.0000E+00 

.0900£4^00 

.000»E^00 


TYPE ICO FOR RESPONSE PLOT 
SINUL I PAUSE 0003 
SinUL SUSP 




figure 1 - 16 (a) 





TUSTIN NUnERATOR COEFriCIEHTS 

-/socE-^ei 

-.l75E+«2 

.25e0E+ei 

-375OE-^02 

.2250E+e2 

.edeeE^oo 

.eooee^de 

.00e®E^0o 

.d00dE-fee 

.0000E+00 

.0O00E^d0 

TtlSTIN DEHOfllHAtOR COEFFICIENTS 

.2812E^01 

-.3aaE>02 

.130SE403 

-.2B2E+03 

.1991E+03 

.0000E+00 

.0000Et00 

.0000E+00 

.0000E+00 

.0000€’^00 

.00e0€'»^00 


TYPE : GO FOR RESPONSE PLOT 
SinUL I PAUSE 0003 
SIftUL SUSP 


Ui 0 



Figure 1-17 (a) 






tustik muher^tor coefficients 

-.i20E>ea 

.TOOOE^dl 

.320dE>d2 

.1750E^e2 

.00O0E4e0 

.0000E'^00 

.0000E>00 

.0000E^00 

.0000E+00 

•eeoeEteo 

TUSTl.'l DENOfllMr rOR COEFFICIENTS 

.5000E+00 

-•600E+01 

.3500E+02 

-•102E+03 

.1125E+03 

.0000E+00 

.C000E+00 

.e000Et00 

.0000Et00 

.0000Ei00 

.0000E+00 


TYPE :GO FOR RESPONSE PLOT 
SINUL I PAUSE 0003 
SinUL SUSP 


f 


Figure 1-18 (a) 



EXACT • SOLID LINE 
TUSTIN - XXX 
SAGE • 4+*^ 


Figure 1-18 (b) 


Figure 1-19 Error Plot Corresponding to Gases 
in Figures 1-14 - 1-18. 
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P/I. SA«P. PER. 


.^a26E-«3 

as47E-e2 

.3304E-^a 

.5362E-ea 

.7403E-ea 


.20OE+0O 

.400E400 

.cedE^ee 

.sooE^ee 

.109Et01 


TVPE iGO FOR ERROR PLOT 
SlflUL t PAUSE 0001 
sinuL SUSP 




Ui 

CTt 


Figure 1-19 (a) 





II. PARALLEL OPERATION OF MICROPROCESSORS FOR REAL-TIME FLIGHT SIMULATION 
2.0 INTRODUCTION 

The use of parallel processing with microprocessors is being studied to 
determine if the new microprocessor technology can yield performance superior 
to that of large digital computers presently being used for real-time flight 
simulation. The simulation cycle time for large digital computers is typically 
about 1/32 second, and this sometimes introduces appreciable phase shift which 
leads to erroneous simulated results for dynamic systems, Hopefully, the use 
of several small processors would allow a decrease in cycle time; but the word 
length would almost certainly be reduced as a practical matter. Thus^ one of 
the fundamental considerations in the study is to determine relationships 
between roundoff error and truncation error, i.e., the relationships between 
small word size with high speed and large word size and slow speed. 

It is assumed, for the present, that high speed can be obtained via par- 
allel processing; but this has not yet been established, The processing meth- 
odology and machine architecture to accomplish this are, themselves, major 
research problems under study. 

As implied above, it is known that the round-off error can become more 
significant when the sampling rate Ci.e,, cycle time, extropolation time, etc.) 
becomes faster. Thus, it is important to determine the best range of sampling 
rates in terms of computational accuracy and stability as a function of trunca- 
tion and round-off error. Once this range of sampling rates is determined, 
ways to implement the simulation on a multi-microprocessor system can be sought 
which will meet a satisfactory sampling rate specification. 

Generally , the truncation error (and the total propagated error, as well) 
in simulation depends not only on the sampling rate and the simulation tech- 
nique used but also on the characteristic frequency of the system to be 
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simulated. Study has been made in the literature [8] on the relative contri- 
bution of the truncation and the round-off error to the error committed at each 
step as a function of system characteristic frequency and the sampling rate, as 
well as the number of bits allotted to the mantissa, A brief review 
follows , 

th 

Let us assume that the j derivative of the state variables of the 
system are given by the relation 


^ = 0 ) 

c 


(2-1) 


where the characteristic frequency o)^ is, in general, complex. A Taylor series 
expansion of x from the point t = nh to t = (n + l)h gives 


1 2 - 1 3 - 

X , - = X + hx + -::r h x + -T” h X + 
n + 1 n n 2 n o n 


• • 


(2-2) 


Substituting (2-1) into (2-2) yields 


12 2 1 3 3 

X - ^ ^ = (1 + ho) + •=• h 0) + T h 0) + • • • )x 

n + 1 c- 2 c 6 e n 


(2-3) 


for which the solution using first-order integration is 


X , -I = (1 + htu )x 
n + 1 c n 


(2-4) 


with the local truncation error 


1 2 2 3 

T, = i h 0) + 0(h ) 

1 Z c 


(2-5) 


The change in x to first order relative to x is 
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(2-6) 


— Ax = ho3 
X n c 

n 


The magnitude of the local relative truncation error is approximately 


" |AxJ 


(2-7) 


where higher-order terms in (2-5) are neglected. Similar expressions can be 
derived for second-order integration, and they are 


— Ax = h03 (1 + ^ h(o ) 
X n c 2 c 

n 


( 2 - 8 ) 


1^2 1 1 


I I ^ I 1 + -7 h03 

n V 2 c 


(2-9) 


The local round-off error R is, when floating-point arithmetic is used, 
bounded by* 


-N, 


-2 



-N. 


X 

= R = 2 ^ 

X 

n 


n 


(2-10) 


where is the number of bits for the mantissa. 

round-off error e is defined as 
r 


The magnitude of the local 


R 


e = 
r 


Ax 


( 2 - 11 ) 


n ' 


*In this case it is assumed that the so-called inherent error can be neg- 
lected with good programming practices. For more detail refer to Ref. [8]. 


60 



Thus, to a first“Order approxiination it has a maximum of 


-N, 


€ = 
r 


(2-12) 


ho3 


and for second-order it has a maximum of 




|ha)^(l + Y ho)^) I 


(2-13) 


It should be noted that the round-off error committed at each step is assumed 
to be random. For purpose of analysis a statistical model is generally adopted 
for the distribution of round-off error. The most common model is to assume 
that the round-off error values are uniformly distributed over the interval 
(0, assumed distribution, the 50 percentile occurs at the mid- 

point, namely (l/2)e^. 

Figure 2-1 shows the local relative truncation and the round-off error 
for first-order and second-order integration expressed by (2-7) and (2-9) , 

(2-12) and (2-13) as functions of hm^. Two values are chosen for one typi- 
cal of a large digital computer and the other a microcomputer with the maximum 
processing capability currently available*. The figure allows a Gonvenient 
comparison of round-off and truncation errors. For example, when the second- 
order Adams-Bashforth (AB-2) method is used on a machine with K - 47 bits and 


^Typical word length of a large digital computer is 60 bits. Out of 60 
bits, 48 bits are assigned for mantissa . Since one bit is reserved for sign of 
the mantissa, the number of significant bits in the mantissa is N-^ = 47. The 
maximum word length of currently available microprocessors in floating-point 
operation is 32 bits. Of the 32 bits the mantissa occupies 24 bits with one 
bit for sign , thus leaving 23 bits in the floating-point representation. 
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y “• 3 

h 03 > 10*" to 10 , the truncation error is at least two orders of magnitude 

c ' " ■ ■ 

~6 “5 

larger than round-off error, while for hw^ > 10 to 10 the situation is just 
the opposite . 

When, at each local step,, one type of error is much larger than the other 
in magnitude, it is considered that the propagated error is dominated by the 
former. Figure 2-1 will be referred to again in a later discussion. 

A practical example has been considered*. The example system is 
described by 

X = Ax + Bu (2-14) 

where 



“0.3236575 

0 

1 ~ 


“0.01785196' 

A = 

0 

0 

1 

, B = 

0 


169521 

0 

-0.4809339. 


^1.379406 _ 


The characteristic roots for the above system are 


= 0 


^2’ ^3 " + jl.07S6 = 1.1512e' 


j (180+69.54)' 


(2-15) 


The step responses using the Euler and second-order Adams-Bashf or th meth- 
ods on a microprocessor with 32-bit floating-point arithmetic have been simu- 
lated by a Hewlett-Packard minicomputer which has 32-bit floating-point arith- 
metic. The system responses obtained have been compared with the solution from 


*This is the short-period approximation for linearized longitudinal air- 
craft dynamic equations of motion, where x^ = 'a (angle of attack), X2 = 0 (pitch 
angle) , and xg = 0 (pitch rate). 
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a fourth-order Runge-Kutta method using h = .001 on a computer with 60-bit 
floating-point arithinetic. This solution is considered to be relatively accu- 
rate. The average absolute relative error is shown in Figs, 2-2, 2-3, and 2-4 
for each state variable, respectively, as a function of sampling frequency. It 
is interesting to note that the total error which is the combinational effect 
of the truncation and round-off error becomes minimum at a certain frequency 
for each state variable. Refer to Fig. 2-1 for an explanation of the results 
in Figs. 2-2, 2-3, and 2-4. 

The magnitude of the characteristic roots for the example system is, from 
(2-15), 

L I = 1.1512 (2-16) 

■ c * 

From Figs, 2-3 and 2-5 the propagated error exhibits a sharp minimum at 

h= 2 msec (2-17) 

Thus, 

Ihm 1 = 0,0023 (2-18) 

From Fig. 2-1, loci 1 and 6 Intersect at hw^ = . 01, i.e. , the trun- 
cation error with the second-order integration method and the maximum routid-off 
error on the 23-bit machine committed at each step would be almost the same in 
magnitude as at - .01, thus competing with each other for contribution to 
the propagated error. Because the slopes of the two loci 1 and 6 are equal in 
magnitude but opposite in sign, we could reasonably have expected that the com- 
binational effect would become minimum around the intersection point of line 1 
and 6- It is noted that the value for in (2-18) is farily close to, but 
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somewhat smaller than, the value for |hm^| at the intersection. This seems to 
be explained, however, by the fact that is the maximum round-off error and 

the round-off error can be any value between zero and at random. Figure 2-1 
also explains why the responses with the AB-2 method on machines with two 
largely different word sizes have almost same accuracy above h - 5 msec. 

^ Particular attention has been paid to the case using a large computer 
with the AB-2 method and h = 1/32 sec, because it was reported [ 8 , 9] that 
this yields real-time flight simulation of satisfactory accuracy except for 
high-performance aircraft with high angular rates. Figures 2-2, 2-3, and 2-4 
show that the average total propagated error with the AB-2 method on a machine 
with 32-bit floating-point arithmetic with a sampling interval h - .002 sec is 
almost two orders of magnitude smaller than that with the AB-2 method on a dig- 
ital machine having 60-bit floating-point arithmetic and h - 1/32 sec. That 
means that ^ if the sampling frequency can be increased to 500 H 3 by means of 
efficient parallel operation of microprocessors, the accuracy of the simulation 
can be significantly improved over the accuracy currently obtained in real-time 
flight simulation on a large digital computer. 

The above discussion was for a single open-loop system. The next step in 
the research will include closed-loop operation and generalization to other 
systems. It has not yet been established that the AB-2 is a suitable method 
for parallel processing; so it must be studied, too. 


67 


REFERENCES 


Parrish, E. A. Jr., Cook, G. , and McVey, E. S., "Real-Time Flight 
Simulation Methodology," Semi-annual Report No. UVA/528085/EE76/103, 

August 1976. 

Fletcher, R. , and Powell, M. J. D., "A Rapidly Convergent Descent Method 
for Minimization," Computer Jnl. , Vol. 6, No. 2 (1963), pp. 163-68. 

Athanassopoulos, J. A., and Warren, A. D., "Design of Discrete Time 
Systems by Mathematical Programming," Proc, 1968 Hawaii Int. Conf. Syst. 
Sci., University of Hawaii Press, Honolulu, 1968, pp, 224-227. 

Burrus, C. S., and Parks, T. W. , "Time Domain Design of Recursive Digital 
Filters," IEEE Trans. Audio Electroacoust. , Vol. AU-18, No, 2, pp. 137-141, 
June 1970. 

Steiglitz, K. , "Computer Aided Design of Recursive Digital Filters," 

IEEE Trans. Audio Electroacoustic , Vol. AU-18, No. 2. (June 1970), 
pp. 123-129. 

Kuo, B. C., "Analysis and Synthesis of Sampled Data Control Systems," 
Prentice-Hall, 1963, 

Beck, M. S., "Adaptive Control— Fundamental Aspects and Their Application," 
Proc. of 1st Anl. Advanced Control Conf. (Dun-Donnelly Publishing, 

Purdue University, April 29-May 1, 1974). 

Wilson, J. W. , and Steinmetz, G. G. , Analysis of Numerical Integration 
Techniques for Real-Time Digital Flight Simulation, NASA TN D-4900, 1968. 

Barker, L. E. , Bowles, R. L., and Williams, L. H. , Development and Applica- 
tion of a Local Linearization Algorithm for the Integration of Quaternion 
Rate Equations in Real-Time Flight Simulation Problems , NASA TN D-7347 , 1973 



